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Abstract 

We propose in this paper a generalized Perron complementation method for uncoupling a consistent linear 
system which involves an irreducible, either singular or nonsingular, M-matrix. We show that this uncoupling 
arises naturally from a regular splitting, which also leads to an efficient iterative scheme for solving the linear 


system. 
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Introduction 


Linear systems involving M-matrices, i.e. whose coefficient matrices are M-matrices, arise in a wide variety of 
fields such as finite Markov chains and partial differential equations. Such linear systems, therefore, have been 
investigated extensively in literature. A majority of the solution methods are iterative, because many applica- 
tions lead to large, often sparse, linear systems; see [5, 6, 9, 11, 22, 24, 28] for Jacobi, Gauss-Seidel, SOR, 
and some preconditioned or accelerated variants, and see [2, 7, 12, 18, 23] for multi-splitting and two-stage 


methods. 


Another useful approach to dealing with large linear systems is divide and conquer methods. In particu- 
lar, Meyer showed in [15] that the Perron eigenproblem can be uncoupled into smaller subproblems via the 
so-called Perron complementation. Such an eigenproblem is indeed a special case of a linear system which 
involves an M-matrix. The Perron complementation technique was later applied to the computation of mean 
first passage times for Markov chains [8] and group generalized inverses of certain Q-matrices [2 1].! However, 


the Perron complementation over linear systems featuring M-matrices has yet to be addressed. 


In this paper, we show how Perron complementation emerges naturally on a linear system involving an 
M-matrix through some regular splitting. This results in an uncoupling of the linear system; furthermore, the 


regular splitting leads to an efficient iterative scheme for solving the linear system. 





1 See [24] for the definition of Q-matrices. 
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Main Results 


For any square matrix X, we shall denote by o (X) and p(X) the spectrum and the spectral radius of X, 


respectively. 


Throughout this paper, 4 = [a;,;] represents an nxn irreducible M-matrix. We mention here a few relevant 
characterizations of such a matrix, which are necessary for the subsequent development. First, there exists an 


irreducible nonnegative matrix B such that 


Az=rI-B, (1) 
where J is an identity matrix and r > p(B); especially, Æ is nonsingular if and only ifr > p(B). Second, 
p(B) > 0 due to the celebrated Perron-Frobenius Theorem (see, for example, [16]). Third, 4 € &, ice. 
aij < 0 for alli # j. Fourth, any (proper) principal submatrix of Æ is a nonsingular M-matrix. Last, the 
inverse of a nonsingular M-matrix is nonnegative. For detailed background material on M-matrices, we refer 
the reader to [1]. 

Consider now a consistent linear system which involves the above Æ in the form 
Ax =b, (2) 
where b € R(A), the range of A. In what follows, we shall focus on a new uncoupling algorithm, along with an 


effective iterative scheme, for solving the linear system (2). 


We assume first that Æ is already expressed as in (1). Partition B into 


, (3) 








where the diagonal blocks are all square, with B,,; being k x k. According to [19], the generalized Perron 
complement of By; in B is defined to be 


G = Boo + Bo, 1 (rl — By,1)'By. (4) 


As shown in [14], G is an irreducible nonnegative matrix with p(G) < p(B), where equality occurs if and 
only ifr = p(B). Thus, F = rI -G is again an irreducible M-matrix, but of a smaller size; moreover, F 
is nonsingular exactly when A is so. The notion of Perron complementation turns out to play a key role in 


uncoupling the linear system (2). 


We point out in passing that the generalized Perron complement can be similarly defined on any principal 
submatrix of B. For this reason, results obtained from the 2 x 2 partitioning of B can be extended to include an 
arbitrary principal submatrix of B or the case when B is partitioned in a general m x m block form with square 


diagonal blocks. 


For any square matrix X, X = M —N, with M being nonsingular, is called a regular splitting when M7! > 0 
and N > 0; see, for example, [27]. Partition Æ in conformity with (3) as 








A= ri-Bii -Bı 9 
7 -Bg, rl- Bog 
Let 
rl — Bı 1 0 0 Bı 2 
M= á dN=M-A= > 5 
-B92 1 rI n 0 Boo ©) 
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It can be easily seen that N > 0 and 


M= >0 


= , 


(rI - By)! 0 
r-'Bo i (rl - By,1)7! roll 


showing that (5) defines a (nontrivial) regular splitting of 4. This regular splitting leads to the following un- 
coupling of the linear system (2). 


by 


THEOREM 1 Let the linear system (2) be partitioned in accordance to (3), with b = | h 
2 





. Suppose that x = | “l | 
x 


is a solution to (2). Then xı and xg satisfy 
Fro =c9, (6) 


where F = rI —G, with G being the generalized Perron complement (4) of B1 ,ı in B, and 
co = B2, 1 (rI — By,1)'b1 + be, 


and 


(rI — By,1)x1 = By 9x9 + b1. (7) 
Proof: By the preceding regular splitting of A, (2) is equivalent to 
(-M7!N)x = M~'b. 
From (5), we further obtain 


I -(rI -B,1)7'Bi2 
0 I-r0 


(rI — By) 1b 
r!Bo (rI — By1)71by +r—!bg 











which yields (6) and (7). 











In the sequel, we shall refer to (6) and (7) as the reduced and the companion linear systems, respectively, 


for (2). Note that (6) is also a linear system involving an irreducible M-matrix. 


In particular, on setting b = 0 and r = p(B), the linear system (2) translates into the (right) Perron eigen- 
problem 


Ba = p(B)z. 
Accordingly, Theorem | reduces to 
COROLLARY 1 [15, Theorem 2.1] Given the irreducible nonnegative matrix B as partitioned in (3) with its Perron 


eigenvector x being conformably partitioned as x = l |, iz follows that 
T2 


Gz = p(B)zxə. 


Theorem | suggests that the solution to (2) can be computed in two separate phases — first on the 
(n — k) x (n — k) reduced linear system (6), and then on the k x k companion linear system (7). However, one 
concern regarding this uncoupling approach comes from its complexity. To address this concern, we estimate as 
follows the amount of work for the case when k ~ n/2, assuming that the solution method is LU factorization. ° 





2 See [4, 10] for LU factorizations of nonsingular or irreducible singular M-matrices. 
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(1) To form the generalized Perron complement G, we first solve 
Y (rI — Bi,1) = Ba,1 


by an LU factorization, which requires 2k? /3 + 2k?(n — k) operations, and we then calculate Y B4 2, which adds 
an extra 2k?(n — k) operations. 

(2) The LU factorization on the reduced linear system (6) needs 2(n — k)?/8 operations. 

(3) The matrix Y and the LU factorization of rI — B1, are stored to be used in the right-hand side of the 
reduced linear system (6) and in the solution of the companion linear system (7), respectively. 


The above two-phase approach, therefore, requires a total of 
t = 2k? /3 + 2(n — k)? /3 + 4k? (n — k) = 2n? /3 — 2k(n — k)(n — Qk) 


operations. In particular, t < 2n°/8 when k < n/2. This explains that, in view of efficiency, we should choose 
a k-value not exceeding n/2. Consequently, the two-phase method can be expected to be less costly than a 
direct LU factorization on the entire linear system. 


The two-phase solution process may be deployed recursively on reduced linear systems. To compute x9 
from (6), for example, we may partition G into 2 x 2 blocks similar to (8), then set up and solve the reduced 
and the companion linear systems for (6). 


We remark that the reduced linear system (6) is also suitable for computing selected entries in the solution 
to the linear system (2). Let œ be a nonempty (proper) subset of (n) = {1, . . . , n}, arranged in ascending order. 


The generalized Perron complement on {n)\q@ in (n) is similarly defined as 
G = Bla, a] + Bla, (n)\a] (rI - B[(n)\q, (n)\a])'B[Xn)\q, a], (8) 


where B[a, £] is the principal submatrix of B on rows œ and columns £. A similar notation is adopted for 


vectors, i.e. x[@] represents the entries of x indexed by a. Then, according to Theorem 1, we arrive at 


COROLLARY 2 Given a, a nonempty subset of (n) = {1,...,n}, and G as in (8), the entries x[a] in the solution 
to (2) can be determined by 


(rI —G)x[a] = bla] + Bla, (n)\a] (rT - B[(n)\a, (n)\a]) "6 [(n)\a]. 


In addition to the uncoupling of the linear system (2), the regular splitting in (5) provides a way of com- 
puting the solution to (2) by iterations. This approach applies to subsequent reduced linear systems as well 
if the two-phase process is used recursively. We proceed next to such an iterative scheme and the issue of its 


convergence. 


The iterative formula on (2) induced naturally via the regular splitting in (5) is given by 
GON Se +c, i=0,1,..., (9) 


where H = M~!N and c = M7!b. It is known that H is convergent, i.e. lim; H’ = 0, if and only if Æ is 
nonsingular [27, Theorem 38.29]. In this case, p(H) < 1; moreover, the asymptotic rate of convergence is 
given by -In p(H). The conclusion below indicates that the iteration (9) converges faster on reduced linear 
systems than on the entire linear system (2). 


THEOREM 2 Suppose that A is nonsingular. Let H4 be the iteration matrix which is induced by the regular splitting 
(5) on the linear system (2). Let Hp be the iteration matrix which is induced by any subsequent regular splitting of the 
same type on the reduced linear system (6). Then 


-1n p(Hp) > -In p(H4). 
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Proof: From (5), we have 


H,= 
4a” | 9 rg 





0 (rl -By1)'Bigs | 


Hence p(H4) =r7!p(G) < 1. 
On the other hand, we partition G into 2 x 2 blocks as 


Gi G12 
G21 Go 








where Gj,; is m xm such that 1 < m < n —k. Using a similar regular splitting as in (5) but on the reduced linear 
system (6), we obtain 
0 (rI-Gi1)'Gi.2 


Hy; = 
Fojo r-IK 








where K is the generalized Perron complement of G1,1 in G. This implies p(Hp) = r7! p(K) < 1. 











Hence the conclusion follows due to the fact p(K) < p(G). 





With the nonsingular case being clarified, we now turn to the scenario when A is singular. 


For a square matrix X, denote 6(X) = max{|A|: 2 € o(X), A # 1}. The index of X, ind(X), is defined to 
be the smallest nonnegative integer m such that rank(X”*!) = rank(X”). The convergence of (9) can then be 


characterized as: 


LEMMA 1 ([1, p. 198]) The iteration (9) converges for any initial x) if 
(1) (H) < land 
(2) ind - A) = 1. 


A matrix H satisfying the two conditions in Lemma | is said to be semiconvergent; besides, lim} sæ H? 
exists. With the regular splitting of 4 in (5), we show next that H = M7!N is indeed semiconvergent when 4 


is singular. 
THEOREM 8 For the splitting M and N as defined in (5), when A is singular, H = M~!N is semiconvergent. 


Proof: According to (5), we have 


0 (rT -By1)'Bis 
H= oe 
| 0 r'G 





where r = p(G). Thus o (H) = a(r7!G) U {0}. The conclusion 6(H) < 1 follows from the irreducibility of 


the generalized Perron complement G. 


On the other hand, 
I -(rI -By1)7'By 


I-H= 
E f=7°C 





By [17, Lemma 1], ind(/ —r7!G) = 1. In addition, from [8, Theorem 7.7.2], 


max{ind(/), ind — r7!G)} < ind(I - H) < ind(J) +ind(I -r7'!G), 
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which yields ind - H) = 1. 


Lemma 3 and Theorem | guarantee that the iteration formula (9) converges to a solution of the linear 
system (2) when A is singular. In the proof of Theorem 1, we observe that the asymptotic rate of convergence 
-In (H) (see [1]) depends on a subdominant eigenvalue of G, i.e. 6(G). Considering a similar situation as in 
Theorem 2, we comment that, in general, 6(K) may not be smaller than 6(G). However, in the special case 
when A is also a Q-matrix with zero row sums, it is shown in [20] that 


Z(K) < Z(G), 


where Z(-) is the coefficient of ergodicity [25] which bounds subdominant eigenvalues; hence the asymptotic 
rate of convergence can be expected not to degenerate much as the solution process is implemented using 
reduced linear systems. 


There is one remaining question which concerns the choice of r and B given the linear system (2), as the 
uncoupling and the iteration processes both hinge on these quantities. To answer this question, we have 
THEOREM 4 Let A = [a;,;], which is an n X n irreducible singular M-matrix. Set r = max; ai i. Then 

(D aii > 0 for alli. 
(2) B =rI — A is irreducible and nonnegative. 
(3) r= p(B) 


Proof: By [1, Theorem 6.4.16], there exists some x > 0 such that Ax = 0. Therefore, for any i, 


iiti + > i jXj = 0. 


j+i 


It is clear that a;,; > 0. If a;,; = 0, then a;,; = 0 for all j # i, which is a contradiction to the irreducibility of A. 


Partition Æ as 








where w = an,n and, without loss of generality, w = max; a; i. Due to the singularity of 4, 





w =v A] u. (10) 
Obviously, B = wI — A > 0 and it is irreducible. It remains to show that w = p(B). Note that 
wl-—A, 1 —u 
B= Í ; 
-yT 0 | 
: : r s i : > s RA P -A7 u 
which implies that w is indeed an eigenvalue of B; in addition, on letting y = 1 , we see that By = wy. 





Since y > 0 but y # 0, it follows, according to [1, Theorem 2.1.4], that y must be an eigenvector associated 
with p(B), i.e. w = p(B). 














When 4 is nonsingular, the conclusions (1) and (2) in Theorem 4, together with r > p(B), follow directly 
from [1, Theorem 6.2.8]. However, the case involving an irreducible singular M-matrix does not seem to be 
stated explicitly in literature except for the case of irreducible Q-matrices.’ Our proof here is different in that 
no Q-matrix argument is needed. 





3In fact, the redundant condition a;,; > 0 appears in some existing work; see [2, p. 309]. 
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Algorithm and Examples 


We begin with a description of the uncoupling algorithm for computing a solution to the linear system (2). 
The situation of a two-level uncoupling is illustrated here. The first level is used to construct the reduced and 
the companion linear systems, while the second level is used to carry out the iterative process for solving the 
reduced linear system. 


Write the linear system (2) as A 2 = 6, with A = [a£ A and bC?) being partitioned conformably as 
AD AWV a) 
by 


AY = L, 
ad) gü 
Ay) Ay 








) 
2 
Choose r = max; al x Set B®) = rI — A, partitioned accordingly as 

ad) a) 
B) 1 Bì 2 


O pO 
Boy By 9 








Next, formulate B®), the generalized Perron complement of B a in BO, 
(2) rD a) d) -1 pQ) 
B`” = By a + By i (rl - By 3) Big 
This leads to the reduced linear system for (2): 
AD =b, (11) 


where 4 =rI - B® and b” eB (rI -B +b. 


Continuing, partition B®? and 6) in conformity as 





D gO (2) 

Bì 1 By 2 (2) _ bi 

p® p® | and B=) say 
2,1 2,2 2 


Similarly, formulate B(?), the generalized Perron complement of B a in B®); 
(8) R2 (2) — R(2))-1 p(2) 
B = Boot Boy (rl Bry) Bio 


The iteration procedure on (11) can then be expressed as 


cD) = HA, +, i=0,1,..., (12) 
where 4 ia 
DnE 
H® 0 (1 - By) Bio 
0 ri B63) 
and 


(2) 





(rI — BO) 
wl 1 
= 2 2)\-1),(2 = 2 
PB (rI - BEP +r- 


Finally, iterate (12) until for some å, |la“*)) — 2 ||/||xf*P || is within a prescribed tolerance, where || - || 
denotes some vector norm. Let x9 = x+). Then a solution to (2) is given by x = : | , with 
T2 





1)\- 1 1 
ay = (rT — BY)! (BY 3a + b}”). 
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(2) 
1,1 


1/8 of the size of A. In a similar fashion, when a three-level uncoupling is used, we choose all BY, to be 


Based on the analysis following Corollary 1, it is advisable to choose both B®? and B,”, to be approximately 


approximately 1/4 of the size of A). 


To demonstrate the above algorithm, we now provide two examples. In both examples, || - Ilo is adopted, 
the tolerance is 10~®, and the initial approximation is set to be a zero vector. The computation is carried out 


in Matlab. 


Example 1. Consider the following n x n irreducible nonsingular M-matrix [9]: 


ie ee ae: om a ip gi 


A= a a3 |? 
a3 a2 
&ı 


ag a} ag ag 1 


with ay = —1/n, ag = -1/(n + 1), and ag = -1/(n + 2). The right-hand side of the linear system (2) is 
generated in such a way that the solution is x = [1,2,..., 2]. We summarize in Table | the results from 
a three-level implementation, which indicate that, in terms of iteration numbers, the three-level algorithm is 
on a par with the adaptive Gauss-Seidel (AGS) method in [26], yet it performs better as compared with the 
Gauss-Seidel method and the modified Gauss-Seidel method in [5] — see a comparison in [9]. 








n number of iterations (three-level) number of iterations (AGS) 
20 37 35 
30 49 50 
50 77 79 
100 145 148 





Table 1: The iteration numbers from a three-level algorithm on Example 1 are reported here. The data from 
the AGS method are quoted from [9]. 


Example 2. This example arises from the discretization of the partial differential equation 


ðu ðu 
Go ane 
Ox? Ay? 





+ 0.50" = f(x,y), (a,y) €Q=(0, 1) x (0, 1), 


with periodic boundary conditions [13]. For this example, the matrix A is given by, in an m x m block form, 


D -I -I 
-I D -I 
A= ; 
-I D -I 
-I -I D 
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where D, with a, = 1+0.5/m and a_ = 1 — 0.5/m, is an m Xx m matrix 


4 -Q4 -a_ 
-a 4 -Ay 
D= 
-a_ 4 -Q4 
—Q4 -a_ 4 


Thus 4 is an m? x m? irreducible singular M-matrix. We choose the right-hand side vector in the linear system 
(2) by b = Ax using a random vector x. The number of iterations from, again, a three-level implementation 
and the residual ||4x* — b|lo, with x* being the computed solution to (2), are summarized in Table 2. 








m number of iterations (three-level) residual 

5 14 5.67 x 1077 
10 55 1.88 x 10-6 
15 110 1.85 x 1078 





Table 2: The iteration numbers and the residuals here are obtained from a three-level algorithm on Example 


2, 
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